Quasiseparable Hessenberg reduction of real diagonal plus low rank matrices and applications

We present a novel algorithm to perform the Hessenberg reduction of an $n\times n$ matrix $A$ of the form $A = D + UV^*$ where $D$ is diagonal with real entries and $U$ and $V$ are $n\times k$ matrices with $k\le n$. The algorithm has a cost of $O(n^2k)$ arithmetic operations and is based on the quasiseparable matrix technology. Applications are shown to solving polynomial eigenvalue problems and some numerical experiments are reported in order to analyze the stability of the approac

Solving polynomial eigenvalue problems by means of the Ehrlich-Aberth method

Given the $n\times n$ matrix polynomial $P(x)=\sum_{i=0}^kP_i x^i$, we consider the associated polynomial eigenvalue problem. This problem, viewed in terms of computing the roots of the scalar polynomial $\det P(x)$, is treated in polynomial form rather than in matrix form by means of the Ehrlich-Aberth iteration. The main computational issues are discussed, namely, the choice of the starting approximations needed to start the Ehrlich-Aberth iteration, the computation of the Newton correction, the halting criterion, and the treatment of eigenvalues at infinity. We arrive at an effective implementation which provides more accurate approximations to the eigenvalues with respect to the methods based on the QZ algorithm. The case of polynomials having special structures, like palindromic, Hamiltonian, symplectic, etc., where the eigenvalues have special symmetries in the complex plane, is considered. A general way to adapt the Ehrlich-Aberth iteration to structured matrix polynomial is introduced. Numerical experiments which confirm the effectiveness of this approach are reported.Comment: Submitted to Linear Algebra App

On Functions of quasi Toeplitz matrices

Let $a(z)=\sum_{i\in\mathbb Z}a_iz^i$ be a complex valued continuous function, defined for $|z|=1$, such that $\sum_{i=-\infty}^{+\infty}|ia_i|<\infty$. Consider the semi-infinite Toeplitz matrix $T(a)=(t_{i,j})_{i,j\in\mathbb Z^+}$ associated with the symbol $a(z)$ such that $t_{i,j}=a_{j-i}$. A quasi-Toeplitz matrix associated with the continuous symbol $a(z)$ is a matrix of the form $A=T(a)+E$ where $E=(e_{i,j})$, $\sum_{i,j\in\mathbb Z^+}|e_{i,j}|<\infty$, and is called a CQT-matrix. Given a function $f(x)$ and a CQT matrix $M$, we provide conditions under which $f(M)$ is well defined and is a CQT matrix. Moreover, we introduce a parametrization of CQT matrices and algorithms for the computation of $f(M)$. We treat the case where $f(x)$ is assigned in terms of power series and the case where $f(x)$ is defined in terms of a Cauchy integral. This analysis is applied also to finite matrices which can be written as the sum of a Toeplitz matrix and of a low rank correction

Efficient cyclic reduction for QBDs with rank structured blocks

We provide effective algorithms for solving block tridiagonal block Toeplitz systems with $m\times m$ quasiseparable blocks, as well as quadratic matrix equations with $m\times m$ quasiseparable coefficients, based on cyclic reduction and on the technology of rank-structured matrices. The algorithms rely on the exponential decay of the singular values of the off-diagonal submatrices generated by cyclic reduction. We provide a formal proof of this decay in the Markovian framework. The results of the numerical experiments that we report confirm a significant speed up over the general algorithms, already starting with the moderately small size $m\approx 10^2$

Shift techniques for Quasi-Birth and Death processes: canonical factorizations and matrix equations

We revisit the shift technique applied to Quasi-Birth and Death (QBD) processes (He, Meini, Rhee, SIAM J. Matrix Anal. Appl., 2001) by bringing the attention to the existence and properties of canonical factorizations. To this regard, we prove new results concerning the solutions of the quadratic matrix equations associated with the QBD. These results find applications to the solution of the Poisson equation for QBDs

Numerical computation of the roots of Mandelbrot polynomials: an experimental analysis

This paper deals with the problem of numerically computing the roots of polynomials $p_k(x)$, $k=1,2,\ldots$, of degree $n=2^k-1$ recursively defined by $p_1(x)=x+1$, $p_k(x)=xp_{k-1}(x)^2+1$. An algorithm based on the Ehrlich-Aberth simultaneous iterations complemented by the Fast Multipole Method, and by the fast search of near neighbors of a set of complex numbers, is provided. The algorithm has a cost of $O(n\log n)$ arithmetic operations per step. A Fortran 95 implementation is given and numerical experiments are carried out. Experimentally, it turns out that the number of iterations needed to arrive at numerical convergence is $O(\log n)$. This allows us to compute the roots of $p_k(x)$ up to degree $n=2^{24}-1$ in about 16 minutes on a laptop with 16 GB RAM, and up to degree $n=2^{28}-1$ in about one hour on a machine with 256 GB RAM. The case of degree $n=2^{30}-1$ would require higher memory and higher precision to separate the roots. With a suitable adaptation of FMM to the limit of 256 GB RAM and by performing the computation in extended precision (i.e. with 10-byte floating point representation) we were able to compute all the roots in about two weeks of CPU time for $n=2^{30}-1$. From the experimental analysis, explicit asymptotic expressions of the real roots of $p_k(x)$ and an explicit expression of $\min_{i\ne j}|\xi_i^{(k)}-\xi_j^{(k)}|$ for the roots $\xi_i^{(k)}$ of $p_k(x)$ are deduced. The approach is extended to classes of polynomials defined by a doubling recurrence